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ABSTRACT 

We explore scenarios for the origin of two different density planets in the Kepler 36 
system in adjacent orbits near the 7:6 mean motion resonance. We find that fine 
tuning is required in the stochastic forcing amplitude, the migration rate and planet 
eccentricities to allow two convergently migrating planets to bypass mean motion 
resonances such as the 4:3, 5:4 and 6:5, and yet allow capture into the 7:6 resonance. 
Stochastic forcing can eject the system from resonance causing a collision between the 
planets, unless the disk causing migration and the stochastic forcing is depleted soon 
after resonance capture. 

We explore a scenario with approximately Mars mass embryos originating exterior 
to the two planets and migrating inwards toward two planets. We find that gravita- 
tional interactions with embryos can nudge the system out of resonances. Numerical 
integrations with about a half dozen embryos can leave the two planets in the 7:6 
resonance. Collisions between planets and embryos have a wide distribution of impact 
angles and velocities ranging from accretionary to disruptive. We find that impacts 
can occur at sufficiently high impact angle and velocity that the mantle of a planet 
could have been stripped, leaving behind a dense core. Some of our integrations show 
the two planets exchanging locations, allowing the outer planet that had experienced 
multiple collisions with embryos to become the innermost planet. A scenario involv- 
ing gravitational interactions and collisions with embryos may account for both the 
proximity of the Kepler 36 planets and their large density contrast. 



1 INTRODUCTION 



Numerical integrations by Deck et al. (2012) illustrate that 



The Kepler mission ( Borucki et al.|2010 l has detected over 
2300 planet candidates (Bat alha et al.|2013[ ) and about one 
third of Keplers candidates are associated with compact 
multiple transiting systems ( Lissauer et al.|2011 l. A statis- 
tical analysis, focused on the probability that binary stars 
are the most likely contaminant, find that most of the mul- 



tiple planet candidates are real planetary systems (Lissauer 
et al. 2012) and that approximately 20% of the total num- 
ber of planet candidate systems are multiple planet systems 
( Latham et al.pOll Lissauer et al.|2011 i. 

Among the newly discovered planetary systems is the 
exotic Kepler 36 system, that hosts two transiting planets 
Kepler-36b and c, with orbital periods of 13.8 and 16.2 days 
and masses of m b = 4.1 and m c = 7.5M®, respectively 
( Carter et al.|2012 i. Of all the multiplanet systems, Kepler- 
36 has the smallest fractional separation between any pair 
of adjacent orbits, and this pair also has one of the largest 
planet density contrasts ( Carter et al.|2012 1 . The estimated 
density for Kep36c is only p c = 0.89 g cm -13 and that for the 
less massive Kep36b is high at pt = 7.46 g cm -3 . The two 
planets' orbits are separated so that they are just exterior to 
the 7:6 first order mean motion resonance. Numerical inte- 



the system is chaotic likely due to the proximity of the 7:6 
resonance and a weak, but high order, two-body resonance. 
Because the Kepler systems are compact, each planet 
fills a sizable fraction of its Hill sphere. This is illustrated 
in Table [I] showing ratios of planetary radii to semi-major 
axes and Hill radii. The distance between the planets semi- 
major axes in units of the outer planet's radius is only 
(a c — a b )/R c = 84.9 and only 4.8 mutual Hill radiiQ As 
explored by Deck et al. (20121, the planets are so close that 



they exhibit chaotic evolution. Because of the planets' prox- 
imity and large size compared to their semi-major axis, it 
is important to check for collisions during numerical inte- 
grations while exploring scenarios for the origin of the sys- 
tem. Because the planets have short orbital periods, their 
orbital velocities are high at approximately 91 and 86 km/s 
for Kep 36b and 36c respectively. The orbital velocities can 
be compared to the escape velocities from the planets that 
are approximately 18.8 and 16.1 km/s. Collisions between 
bodies in this system would be at high velocity compared 
to the planets' escape velocities. A scenario involving a high 
velocity collision, such as proposed for the formation of Mer- 



grations by Carter et al. (20121, constrained by transit tim- 



ing observations, not only have measured the planet masses, 
but restrict the planet orbital eccentricities to less than 0.04. 



1 The mutual Hill radius r m n = 



' m b +m c 
v 3M, 



)(-4-*). 
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Table l. Some dimcnsionless quantities in the Kepler 36 
planetary system 



Planet mass/Stellar mass 
Orbital velocity/Escape velocity 
Semi-major axis/Hill radius 
Hill radius/Planet radius 
Semi-major Axis/Planet radius 

These quantities arc computed using masses, planetary radii and 
orbital periods by|Carter et al.|l|2012[). 



Planet b 


Planet c 


.15 X 1CT 5 


2.09 X 10~ 5 


4.8 


5.3 


63.9 


52.3 


29.0 


16.0 


1852 


838 



cury, and accounting for Mercury's high density ( Benz et al 



2007), could account for the high density of Kep 36b. 



In this paper we investigate scenarios for the origin of 
the Kepler 36 system, focusing not only on the proximity of 
the two planets but on the possibility of collisions between 
bodies. 



2 FINE TUNING IN THE STOCHASTIC 
MIGRATION SCENARIO 

Previous works have explored the properties of planet semi- 
major axis distributions that are predicted from stochas- 



tic migration models (Rein 20131. It is assumed that the 



semi-major axes and eccentricities of the planets can vary 
to due to torques transferred between embryos or planets 
to a gaseous disk via driving of density waves ( Ward|[l997 



Tanaka et al.|2002 l. The gas disk is no longer present, hav- 
ing long since dissipated. A two planet system near the 7:6 
was not predicted by the stochastic migration scenario and 
even high mass planet pairs in the 4:3 resonance are diffi- 
cult to explain (Rei n et al.|2012 |. We first explore whether a 
stochastically migrating two planet system would be likely to 
be subsequently found near the 7:6 mean motion resonance. 
Concurrent to our study is the recent numerical exploration 



by Paardekooper et al. ( 2013 1 who have shown that stochas- 



tic migration is a viable way to allow capture into the 7:6 
resonance for the Kepler 36 planets. Our approach in this 
section is primarily analytical rather than focused on hy- 
drodynamical simulations and so is complimentary to their 
study. 



2.1 Sensitivity to migration rate and initial planet 
eccentricity 

We consider a setting where the two planets approach one 
another due to convergent migration. As the system drifts it 
encounters mean motion resonances between the two plan- 
ets. A first order j : j — 1 resonance is a commensurability 
where jn c w (j — l)nb where nb,n c are the mean motions 
of the two planets. Capture into a mean motion resonance 
is likely if the planet eccentricities are low, the planets ap- 
proach each other, and the relative drift rate is smooth and 
slow (e.g.,|Murray fc Dermott| 1999 Hcnrard 1982; Bor deries] 
fc Goldreich|1984||Quillen||2006||Mustill fc Wyatt||2011[ ). 



Under what conditions can a smoothly migrating sys- 
tem bypass lower j first order mean motion resonances such 
as the 3:2 and 4:3 resonances, but capture in a higher j one 



such as the 7:6 resonance? Higher j first order mean mo- 
tion resonances are stronger, wider and have faster libration 
frequencies than lower j ones. If the migration rate is suf- 
ficiently high, then first order resonances such as the 3:2, 
4:3, 5:4 and 6:5 would be bypassed, but capture into the 7:6 
resonance would still occur. 

Capture of a low mass particle into a first order mean 
motion resonance with a planet can be modeled with a one 



dimensional Hamiltonian system (e.g., Murray & Dermott 



1999 



1990 



Henrard|[l982| |Borderies fc Goldreich|1984| |Malhotra 



Quillen|2006| |Mustill fc Wyatt|2011 i 



H(T, 4>) = AT' + bT + eT l,i cos 



(1) 



where T = ^/a(l — y/1 — e 2 ) is a Poincare momentum for a 
particle with semi-major axis a and eccentricity e. We con- 
sider capture of a massless particle into resonance with a 
planet. We neglect the fact that the Kepler 36 system con- 
sists of two planets with similar masses. To order of mag- 
nitude, this is justified as the more general formulation by 



Mustill & Wyatt ( 2011 1 only differs by factors of order unity 
from the case of a massless particle near a planet. Above 
we have assumed units GM* — 1 with G the gravitational 
constant and M, the mass of the host star. The resonant 
argument, 



: j\ - (j - l)Ap + ZU, 



(2) 



corresponds to a slowly varying angle in proximity to the j : 
j — 1 resonance for an object external to a planet. Here A, A p 
are the mean longitudes of particle and planet, respectively, 
and w is the particle's longitude of pericenter. 

The coefficients, A, b, e, depend upon the integer j and 
the ratio of particle to planet semi-major axes, a = a/a p . For 
a drifting system, the proximity to resonance, described by 
frequency b, is a function of time. This model is derived by 
expanding the Keplerian Hamiltonian near resonance and 
adding the lowest order eccentricity term in the disturb- 
ing function. The coefficients A w — j 2 and e ss Mp/27( a ) 



[Murray fc Dermott 1999 Quillen 2006 Mustill fc Wy- 



att||2011 1 where /27(a) is a function of Laplace coefficients 



and is given in the appendix by Murray & Dermott ( 1999 1 



Here \i v is the ratio of planet to stellar mass. Because the 
7:6 mean motion resonance has a moderately high j value, 
we can estimate the dimensions of the resonance using the 
asymptotic limit; a — > 1. In this limit, the Laplace coeffi- 



M) 



cient, b~ L(pi) ~ ln(l — a) (Quillen 



'1/2 



-j for a 



2011 1, and the coefficient 
0' - !)/j s et by the 



/27(a) ~ -(1 - a)- 1 
resonance condition. 

The particle can only be captured into resonance if the 
drift rate, b, is slower than a critical rate that can be es- 
timated through dimensional arguments and delineates the 
adiabatic limit (|Quillen|[200"6| |Mustill fc Wyatt||2011[>. Be- 



cause the frequency b ~ jn— (j — l)n p is small in resonance, 
the Hamiltonian contains a single dimensional unit of time 



tlib 



-2/3 



A\ 



-1/3 



(3) 



that is also approximately the inverse of the libration fre- 
quency in resonance. The Hamiltonian contains a single di- 
mensional unit of momentum 



2/3 



(4) 
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that is approximately the size scale of eccentricity oscilla- 
tions for a system in resonance at low eccentricity. 

Capture into resonance is only possible when the drift 
rate is below a critical value that is approximately the square 
of the inverse of the libration timescale; \b\ 



;4/3| A |2/3 



corresponding to a critical planetary migration rate that is 
to order of magnitude 



1 r lcrit 



4/3 -5/3 

Up 3 , 



(•») 



where we have used the high j (or a — ¥ 1) limit for the co- 
efficient e. Taking the square root of the momentum scale 
we can estimate the critical eccentricity (below which cap- 
ture has probability 1 in the adiabatic limit; [Bordcrics & 
Goldreich|1984l |Malhotra|1990| l of 



&crit 



1/3 --1/3 

Mp 3 ■ 



(6) 



This eccentricity value also corresponds to a mean value 
for the size of the eccentricity jump that occurs when the 
system crosses the resonance instead of capturing into it 
| |Quillen|200"o]|Mustill fc Wyatt|20l"T| . By differentiating the 
above two equations with respect to j, we can estimate the 
difference between critical drift rates for neighboring first 
order resonances, 



(7) 



1 


CLTlcrit 


5 


Tlcrit 


dj 


3j 


and critical 


eccentricity 


1 


£*>&crit 


1 


&crit 


dj 


3j 



(8) 



For j = 7 the fractional difference in critical drift rates be- 
tween the 7:6 and 6:5 resonances, is of order 1/4 and is 

small. This is consistent with Figure 11 by Mustill 



att 



(2011 1 illustrating critical drift rates for first and second 



order resonances. 

The fraction difference between the critical eccentricity 
allowing capture into the 7:6 but not 6:5 is even more severe; 
it is only 5%. It would be difficult to maintain the system 
at a particular initial eccentricity value as the two planets 
undergo secular oscillations as they migrate inwards, eccen- 
tricity damping is likely during migration and as the sys- 
tem crosses resonances, planets can increase in eccentricity. 
First order mean motion resonances also contain a corota- 
tion term that depends on the argument 



: jX - (j - l)X p + ' 



(9) 



and when the planet's eccentricity is similar to or above 
the critical eccentricity, this term can prevent capture into 



resonance (Quillen 20061. For the Kepler 36 planets, with 



planet to stellar mass ratio listed in Table 111 the critical ec- 
centricity is approximately e cr ,: t ~ 0.03. If the eccentricity 
damping rate is low, then we expect the planet eccentrici- 
ties to increase if the system crosses first order resonances 
without capturing. When the particle eccentricity is initially 
above the critical eccentricity value (r > r res ), the proba- 
bility of capture into resonance is not zero but is reduced. 
Studies of stability have shown that higher eccentricity two 
body systems are often less stable than low eccentricity sys- 
tems (e.g., |Kley et al.|2004| [Mustill fc Wyatt|2011[ ). The 7:6 
resonance is so close to the region of resonance overlap at 
low eccentricity that even moderate planetary eccentricities 
put the system within the chaotic zone. Given the number of 



processes affecting planet eccentricity, it would be extremely 
unlikely that a body could remain near the critical eccentric- 
ity value for the 7:6 resonance as the system drifted inward, 
approaching the resonance. 

We conclude that fine tuning is required to adjust either 
the migration rate or planet eccentricities so that capture 
(by smooth convergent migration) into the 7:6 mean motion 
resonance is possible but not other nearby mean motion res- 
onances. 



2.2 Stochastic migration 

Stochastic forcing of a planet by disk turbulence or planetes- 
imal scattering can prevent resonance capture or kick a body 
out of resonance ( [Zhou et al.||2002| |Murray-Clay fc Chiang 
|Rein fc Papaloizou||2009[ |Ketchum et al.||2011| |Rein 



2006 



2013 



Paardekooper et al.|2013[). Stochastic forcing has been 



described in terms of the variance of the angular momentum 
change per orbit or equivalently in terms of a diffusion coef- 
ficient in angular momentum or semi-major axis. A planet 
wanders in semi-major axis a typical distance 



after a time t, where stochastic forcing is described in terms 
of a diffusion coefficient coefficient, D a . If the diffusion 
causes a random walk in semi-major axis, the above equation 
gives the standard deviation of the distribution of distances 
travelled. For a in units of the planet's semi-major axis, a p , 
and t in units of n p 1 (the inverse of the planet's mean mo- 
tion), D a is in units of a p n p . The diffusion coefficient, D a , is 
approximately equivalent to the square of the a parameter 



adopted by Rein (20131, describing the ratio of the stochas- 



tic force perturbation per orbit in units of the gravitational 
force from the central star. 

A stochastically forced body escapes resonance when 
the particle's semi-major axis varies by an amount of order 
the resonance width ( |Murray-Clay fc Ch iang 2006L Equiva- 
lently, a variation in the coefficient 6 (from our Hamiltonian; 
equation nj greater than the resonance libration frequency 
would let the particle escape resonance; 



Sb~ jSn>t ti t, 



(10) 



for 5n the size of the variation in the mean motion of the 
planet. As long as the planet is at low eccentricity, this im- 
plies that the particle would escape resonance after a typical 
timescale, t esc , with 

jVDj--;.~\e\ 2/3 \A\ 1/3 ~nl /3 f /3 (11) 

giving a timescale for escaping the resonance 
t esc ~ J D-Vp /3 / /3 , (12) 

where we have used equation [3] for the libration timescale. 

2.3 Timescale to escape resonance via stochastic 
forcing when at the balancing eccentricity in 
resonance 

Once captured in resonance, the body's eccentricity in- 
creases with growth rate dependent on the drift rate. The 
body's eccentricity increases until it reaches a value etaiance, 
where there is a balance between eccentricity growth due 
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to the resonance and eccentricity damping from dissipation 



forces (Gomes 1998 Lee & Pcalc 20021. As the resonant 



width increases with eccentricity, the above relation (equa- 
tion 12 1 for the escape timescale is a lower limit as it was es- 



timated at low eccentricity. To estimate the escape timescale 
after capture into resonance we first estimate the balancing 
eccentricity, etaiance, and then we revise our estimate for the 
escape timescale. 

Once captured into resonance, the Poincare momentum 
r increases; equivalently the particle's eccentricity increases. 
Hamilton's equation (using equation nl) gives 



f-=Mr + 6 + £r 



V2, 



In resonance, the resonant argument librates about a fixed 
value so the average {(j>) (averaged over a libration or oscilla- 
tion timescale) is constant and (</>} = 0. As F increases, the 
first two terms on the right hand side of the above equation 
dominate over the third. Using these two approximations, 
the time derivative of Hamilton's equation gives 2AT + b = 
and so an increase rate in the Poincare momentum 



b 
2A' 



(13) 



Using A ~ — j 2 and a low eccentricity approximation for the 
Poincare momentum, the eccentricity growth rate in reso- 
nance is of order 



1 



(14) 



where we have written the migration rate (setting b) in terms 
of a migration timescale 



T a = a/a 
Equation 



(15) 



14 



implies that e 2 oc t in resonance and is consis- 
tent with previous work ( Mustill fc Wyatt|2011 their equa- 
tion 6). Eccentricity damping on a timescale, 



= e/e, 



(16) 



can balance the eccentricity growth rate given by equation 
|14| This occurs at an eccentricity that is estimated by bal- 
ancing the two rates for e; 



G-balanc 



JTa 



where we have used a parameter 

K = T a /T e 



(17) 



(18) 



that describes the ratio of the eccentricity damping to mi- 
gration rate. As is true in the case of dust particles captured 
into resonance by Poynting Robertson drag, for the high j 
resonances, the balancing eccentricity is lower than for lower 
j resonances ( Liou fc Zook|1997 k For high rates of eccentric- 
ity damping (high K), the resonant system does not reach as 
high eccentricities and so is more stable ( |Lee fe Peale 2002 
Ketchumet~aL|201U 

When the limiting or balancing eccentricity in reso- 
nance is higher than that set by dimensional analysis, or 



r >r r 



the Hamiltonian can be approximated bjrl 



H = Ap 2 + b'p + eT ' cos <j> 



where To = (r) is the average value of F, we have defined 
a new momentum, p = V — To, and we have incorporated 
the shift in the distance to resonance with a new coefficient 
b' . The resonant width is set by the mean eccentricity or 
To ~ eLiance/2- In this limit the inverse of the resonant 
libration frequency is 



t ub = l/d\Ae\T, 



■a/2 



.-3/2 -1/2 -1/2 
' J ^P e balani 



(19) 



This gives a timescale to escape resonance (from 5b — 

jVDatesc - t' 1 ) Or 



tesc — 



AeT, 



1/2 



] 2 D a 



J fJ'pG-balance 

D~ a • 



(20) 



Inserting equation W7\ for the balancing or limiting eccentric- 
ity in resonance we estimate an escape timescale 



i 1/2 M P 



This timescale should be a better estimate than the 
timescale given in equation [12] as it takes into account the 
eccentricity of the body in resonance. 

Ignoring the weak dependence on K and j in equation 
|21[ the timescale to escape resonance diffusively is approxi- 
m ately t e sc ~ / x^D" 1 . The stochasticity parameter assumed 
by|Rein||2013| of a 



10 , corresponding to D a ~ 10~ , 
gives an escape timescale only of order ~ 10 6 orbital periods 
or ~ 10 years taking into account the orbital period of the 
Kepler 36 planets. For the stochastic migration to account 
for the proximity of the Kepler 36 planets, the stochastic 
forcing would have had to dissipate on a short timescale, 
otherwise stochastic forcing would have removed the sys- 



tem from resonance. Kley et al. ( 2004 \ pointed out a similar 



fine tuning problem. A system that continues to migrate 
after resonance capture can become unstable as the plane- 
tary eccentricities increase. Here the system that continues 
to be stochastically forced after capture into resonance can 
become unstable. In the situation discussed by |Kley et al.| 
(2004), the gaseous disk, responsible for planetary migra- 
tion must dissipate on a timescale short compared to the 
evolution timescale. Here, the source of the stochastic forc- 
ing must dissipate on a timescale short compared to the 
evolution timescale of the system. 

The escape timescale estimated in equation [21] repre- 
sents an estimate for the time that a stochastically forced 
pair of planets in resonance would escape the resonance. 
This estimate neglects the possibility that dissipative forces 
(such as eccentricity damping) could continually nudge the 
system into a stable region and the possibility that the 
stochastic variations do not cause a true random walk in the 
semi-major axis. These factors might increase the lifetime 
of a resonant but stochastically forced system. Systems in 
resonance may be unstable under gravitational forces alone 
and this might decrease the lifetime of the resonant system. 
We will discuss these issues further after we illustrate the 
behavior of stochastically forced systems numerically. 



2 See the appendix on the Andoyer Hamiltonian by |Ferraz-Mello| 
(2007). 
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2.4 Numerical Integrations 

We explore stochastic migration scenarios by integrating a 
few body system under the influence of gravity (a few planets 
and the central star) and including a Stokes drag-like form 
for dissipation that induces both migration and eccentricity 
damping. The drag gives a force per unit mass in the form 



adopted by Beauge et al. ( 2006 \ 



drag — 



V 

2^~ 



(22) 



where v is the planet velocity and v c is the velocity of a 
planet in a circular orbit at the current radius (from the 
star) of the planet. We use a 4th order adaptive step-size 
Hermite integrator (that described by Makino & Aarseth 
1992 ) with the addition of the above drag force. 



We also introduce random velocity variations in the or- 
bital plane (as did iKetchum et al.|2011 1. The distribution of 
velocity kicks is described by a Gaussian distribution with 
a standard deviation, e^ick, in units of the speed of a parti- 
cle in a circular orbit at the current position of the planet. 
The velocity kicks are given twice per orbit. The diffusion 
coefficient discussed above; D a ~ e\ ick . Migration rates, ec- 
centricity damping timescales and stochastic forcing param- 
eters T a ,T e ,ekick are set individually for each massive body 
in the integration except the central star. 

During the integration, we continually check for col- 
lisions between bodies (planets, embryos or the star). For 
most of our integrations we use the measured masses and 
radii of the planets Kepler 36b and c and its host star. A 
collision is identified when the two bodies have distance be- 
tween their center of masses that is equal to or within the 



sum of their radii. Following Asphaug (20101, at the mo- 



ment of impact, we define the angle of impact as the angle 
between between the relative velocity vector and the vector 
between the two center of masses; 



^impact ' Air 

|v im p act ||Ar| 



(23) 



Here 9% m — for a direct normal collision and dim = 90° 
for a grazing collision where the surfaces barely touch. A 
grazing impact is defined as one with angle such that the 
center of mass of the smaller body would graze the surface 
of the larger body, 

^sin-^), (24) 

where R is the radius of the larger body and r is the radius 
of the smaller body. If 9 im > g then the surface of the larger 
body can be removed during the collision ( Asphaug|2 010). 
We work in units of the innermost planet's initial semi- 
major axis, and the mass of the central star. Time is such 
that the innermost planet's initial orbital period has a value 
of 27r. We primarily plot figures in units of the innermost 
planet's initial orbital period. 



axis is 1.42 times that of the inner planet. The forced migra- 
tion timescale r a = — 10 6 for the outer planet and is infinity 
for the inner planet. For the outer planet K = 20, setting 
the eccentricity damping timescale. Estimates of the ratio 
of the eccentricity damping ratio, K = T e /r e , from hydro- 
dynamic simulations range from K 



1 (|Kley et al 
10 ( jCresswell fc Nelson|2008l ) to 100 (Bitsch fc Kley 



2004), 



2010) 



We have chosen values for K in the middle of this range. 
The stochasticity parameter ekick — 10 -5 and 1.5 x 10~ 5 for 
Figures [I] and [21 respectively. Initial orbital angles were ran- 
domly chosen. The initial inclinations were extremely small 
but nonzero; the outermost planet has an initial inclination 
of 0.1°. 

In Figure [I] and [2] the top panel shows semi-major axes, 
with error bars set according to the planet's eccentricity so 
that they illustrate pericenters and epicenters of the orbit 
as a function of time. Both planets are shown together in 
this panel. The middle panel shows the ratio of orbital pe- 
riods and major mean motion resonances are shown with 
grey horizontal lines. The bottom panel shows the difference 
between the planets longitudes of pericenter and illustrates 
when the two bodies are in resonancerlln both integrations, 
the time between reaching the vicinity of the (or capture 
into) 7:6 resonance and a planet/planet collision was short 
compared to the migration timescale. The two integrations 
shown were chosen from approximately 30 integrations with 
different initial orbital angles planet eccentricities, r a , K and 
e-kick values. The two integrations shown were chosen be- 
cause the two planets remained the longest time near or in 
the 7:6 resonance. When we adjust the stochasticity param- 
eter so that lower j mean motion resonances are bypassed, 
we find that the system is unstable and collisions between 
the two planets are likely, on a short timescale, independent 
of whether the planets are begun at low or high initial ec- 
centricity. The short lifetime of stochastically forced planets 
in the 7:6 resonance we see in the integrations is consistent 
with the short timescale we have estimated analytically for 
the resonant lifetime. 

Both our order of magnitude calculations and our nu- 
merical integrations suggest that capture into the 7:6 reso- 
nance is possible for stochastically forced migrating planets, 
but that the planets are unlikely to remain there for long. 
This is not confirmed by the hydrodynamic study (and ac- 
companying thorough N-body study) by Paardckooper et 
al. (20131 who find that the planets are likely to be placed 



into and remain in stable resonant configurations despite the 
stochastic forcing. It may be possible to reconcile our order 
of magnitude estimates with their results through a more 
careful study of the stochastic process and so improve on 
our simplistic prescriptions of the diffusive process. 



2.5 Illustrations of stochastic migration for two 
planets 

Example stochastic migration integrations for two planets 
are shown in Figures [I] and [2] For both integrations, the 
two planet masses and radii are equivalent to those mea- 
sured for Kep b and c. The outer planet's initial semi-major 



3 The perturbation terms in the Hamiltonain with resonant ar- 
guments cos <f> and cos <j>' have opposite sign. When the different 
between the planets' longitudes of pericenter is W\, — vo c fa it, 
the two terms constructively add and the resonance is effectively 
stronger. 
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Figure 1. Integrations of two planets, where the outer planet 
migrates stochastically with migration timescale r a = 10 6 , ec- 
centricity damping factor K = 20, and stochasticity parameter 
e kick = 10 -5 . The planets have masses and radii of those mea- 
sured for Kepler 36b and c. The inner planet was initially started 
with a moderate eccentricity of 0.2. The top panel shows semi- 
major axes, with error bars set according to the planet's eccentric- 
ity so that they illustrate pericenters and epicenters of the orbit 
as a function of time. The middle panel shows the ratio of or- 
bital periods and major mean motion resonances are shown with 
horizontal lines. The bottom panel shows the difference between 
the planet's longitudes of pericenter. The eccentricity variations 
are due to secular oscillations. The capture probability is reduced 
both by the eccentricity and the stochasticity. If the eccentricity, 
migration rate and stochasticity of the planets is sufficiently high 
then low j resonances can be bypassed. In this simulation the 
planets were captured into the 7:6 resonance, however the system 
did not remain there for long. 



3 STOCHASTIC FORCING BY PLANETARY 
EMBRYOS 

We consider the possibility that stochastic forcing is asso- 
ciated with encounters with planetary embryos, bodies of 
order Mars or Earth mass, rather than turbulence from a 
gaseous disk. A disk edge can trap embryos forming a 'planet 
trap' ( Masset et al. 2006 ) where planetary embryos can col- 
lide with each other ( Morbidelli et al!p 008 ) or with planets 
in the vicinity of the disk edge. A feature of this scenario 
is that collisions between embryos and planets are poten- 
tially capable of stripping the outer layer of a planet and so 
accounting for the high density of Kepler 36b. 

We consider a group of planetary embryos, with planet 
to stellar mass ratio fJ. em bryo = 7x 10~ 7 (equivalent to 2.5 the 
mass of Mars), initially exterior to two planets. This embryo 
mass size scale was chosen so that collisions can be highly 
disruptive (Benz et al. 2007 Asphaug 20101. We integrate 



Figure 2. Similar to Figure [T] except both planets had ini- 
tially low eccentricities of 0.02 and the stochasticity parameter 
e kick = 1-5 x 10~ 5 . After 40,000 orbital periods, the two planets 
experienced an era of strong close encounters that ended with a 
collision. 



turbulence in a gas disk) is set to zero for all bodies. As did 



Morbidelli et al. (20081, we allow the embryos to migrate 



faster than the planets. This was motivated by comparing 
a type I migration rate (appropriate for embryos embedded 
in a disk) to a type II migration rate (appropriate for gap 
opening planets) , within the context of torques transferred 
between embryos or planets to a gaseous disk via driving 



of density waves (Ward 1997 Tanaka et al. 20021. If and 



the system allowing the planetary embryos and the outer- 
most planet to migrate inwards. Stochastic forcing, (due to 



when the two planets exchange locations so the inner planet 
becomes the outer one, then the migration and eccentricity 
damping rates are also swapped. This would be expected 
if the outer planet opened a gap in the gaseous disk, and 
the gas did not significantly penetrate within its orbit. This 
procedure allows convergent migration at all times for the 
two planets. In the event of a collision, momentum and mass 
were conserved and a new larger body created at the loca- 
tion of the center of mass of the collision with a density 
equivalent to the mean density of both progenitors. 

We ran two sets of 10 numerical integrations that we 
denote the X-series and the Y-series. The X-series set has 
planet masses and radii equivalent to those of the Kepler 36 
b and c. The Y-series set has two equal planet masses with 
mass and radius equivalent to Kepler 36c. Parameters for 
these integrations are listed and described in Table [2] 

From the 10 X-series of integrations the end states 
were 4 integrations with collisions between the two planets, 
though at the end of one of these an embryo was trapped 
in the 5:4 resonance with the remaining planet, one integra- 
tion each with planets in the 4:3, 5:4, 7:5 and 6:5 resonances, 
one simulation with the planets near the 7:5 resonance but 
at moderate eccentricity so likely to experience a collision 
later on, and a resonant chain with the two planets in the 
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Table 2. Integrations with Two planets and Embryos 



Integration series 


X-series Y-series 


M6 


1.15 X 10" 5 2.09 X 1CT 5 


He 


2.09 x 10" 5 2.09 x 10~ 5 


Rb 


0.00054 0.001 


Re 


0.001 


[^embryo 


7 x 10~ 7 


^embryo 


0.0003 


T a Outer planet 


2 x 10 6 


r a Embryos 


10 6 


K Outer planet 


20 


K Embryos 


10 


Number of Embryos 


7 



Here fn,, fi c , ^, em hryo refer to the body to stellar mass ratio for the 
inner planet, outer planet and embryo masses, respectively. The 
parameters R^,R C , R^mbryo are the body radii in units of the ini- 
tial semi-major axis of the innermost planet. The embryo radius 
was chosen so that embryos have the same mean density as Ke- 
pler 36c. The innermost planet has no forced migration. There is 
no stochastic forcing in these integrations. The initial semi-major 
axis ratio of the two planets is a c /ai, = 1.4. The embryos were 
initially separated so that the ratio of the semi-axis to that of 
the next nearest body is a^i/a^ = 1.1. Embryos were begun on 
circular orbits exterior to the planets. Initial eccentricities of all 
bodies were set to zero. The parameter r a describes a timescalc for 
migration and is given in units of the inverse of the inner planet's 
initial mean motion. The parameter K is the ratio of the eccen- 
tricity damping rate to the migration rate. Initial orbital angles 
were chosen randomly. Initial inclinations were chosen randomly 
and within 0.1 degrees of the mid plane. Integrations were run 
for a time of between 200000 and 400000 orbital periods (of the 
inner planet initially) allowing all inward migrating embryos to 
interact with the planets. Three of these integrations are shown 
in Figures [3] [2] and [5] 



4:3 resonance and an embryo in a 3:2 resonance with the 
outer planet. In two integrations of these integrations, the 
two planets swapped locations. From the 10 Y-series of inte- 
grations the end states were 3 integrations with two planets 
in the 4:3 resonance, one integration with planets in the 5:4, 
three in the 6:5, one in the 7:5 and one in the 7:6 resonance. 
The final result of one integration was a resonant chain with 
the inner two planets in a 7:5 resonance and the remaining 
embryo in a 4:3 with the outer planet. In the Y-series there 
were no planet/planet collision nor did the planets swap lo- 
cations. 

Figures [3] and [4] show integrations from the X-series and 
Figure [5] shows an integration from the Y-series. The top 
panel shows semi-major axes, with error bars set according 
to the body's eccentricity so that they illustrate pericenters 
and epicenters of the orbit as a function of time. Here the 
red and green points are the two planets and the other color 
points correspond to the planetary embryos. The middle 
panel shows the ratio of orbital periods for the two planets 
(outer divided by inner) and major mean motion resonances 
are shown with horizontal lines. The bottom panel shows in- 
clinations of the bodies, with the same color points as the top 
panel. Inclinations were measured with respect to the initial 
orbital mid plane. These integrations illustrate that encoun- 



ters with planetary embryos can knock the convergently mi- 
grating planets out of resonances and so allow them to get 
closer together. Gravitational encounters with embryos can 
serve as a type of stochastic forcing that pushes the sys- 
tem past resonances such as the 4:3 and 5:4 resonances, but 
sometimes lets the planets into the 7:6 resonance. 

What mass embryo is sufficiently massive to push a 
planet out of resonance? We estimate the size of a varia- 
tion in semi-major axis or equivalently energy that would 
remove a system from resonance. As we did previously, we 
estimate the width of the resonance as 6 ~ jSn < tj ih for a 
variation 5n from the center of the resonance. A variation in 
mean motion can be related to a change in semi-major axis 
and so energy. Using equation [19] for the libration frequency 
(describing the resonant width at an eccentricity et, a ; ance ), 
we estimate a variation in orbital energy of order 

— > i 1/ V / V /2 (25 s ! 

JJ ~J Mp e balance' \ Z0 ) 

would remove a planet from the j : j — 1 resonance with an- 
other planet. Using the impulse approximation we estimate 
that an embryo undergoing a gravitational encounter with 
impact parameter b en (in units of the planet's semi-major 
axis) with a planet would cause a change in orbital energy 
of order 

Oly ^embryo /Oft\ 

1-J VenOen 

where the encounter velocity v en is units of the planet's cir- 
cular velocity. Equating these two relations we estimate that 
an encounter with v en b en and embryo mass ratio ^embryo 



^embryo 



1/2 1/2 1/2 
Mp e bala 



„VenOen 



(27) 



can remove the planet from resonance. The minimum impact 
parameter is given by the radius of the planet and is of 
order 1/1000 in units of the semi-major axis (see Table [l]). 
Taking v en ~ 0.2 the circular velocity (measured from the 
distribution of collisions between planets and embryos as we 
will discuss below), b en = 1/1000, and the mass ratio of Kep 
36c, we estimate that an embryo with mass ratio of order 
(i ~ 10~ 6 is required to knock the Kepler 36 system out 
of a low j first order mean motion resonance. The embryo 
mass chosen for our integrations is consistent with this mass 
estimate. Convergent migration of the two planets allows 
them to be captured in resonances such as the 3:2 resonance 
(see Figure[3|. In most but not all cases, the planets escaped 
resonance because of a collision with an embryo. In some 
cases, close encounters were sufficient to kick the system 
out of resonance. 

We find that there is a diversity of possible end states 
exhibited by the integrations. These include resonant config- 
urations with two planets alone, including high j resonances 
such as the 7:6 (see Figure [El. This end state is in a fairly 
stable region, though longer integrations that include a de- 
pletion timescale to remove the dissipative forces are needed 
to determine if this end state would be long lasting and 
leave the two planets in a configuration consistent with that 
observed for the Kepler 36 system. The Kepler 36 system 
lies just outside the 7:6 resonance, so subsequent evolution 
would have to account not only for the stability of the system 
but removal from the resonance. A couple of our integrations 
ended with a last embryo in a resonant chain configuration 
where each pair of bodies is in resonance (see Figure HI. The 
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Kepler planetary systems include resonant chains (e.g., the 
KOI 730 system Lissauer et al.|2011 I , so the diversity of our 
integrations is not inconsistent with the diversity of Kepler 
planetary systems. Two of the integrations show the two 
planets exchanging locations (see Figure p3j). This is reminis- 
cent of the exchange of Neptune and Uranus in some 'Nice' 



model simulations of the early Solar system evolution (Mor- 



bidelli 20101. Among our 20 integrations, we find that the 



outermost planet experiences approximately twice as many 
collisions with embryos as the innermost planet. These colli- 
sions could have stripped the envelope of this planet, leaving 
it at higher mean density than the other one. Because the 
two planets swap locations, the end result might be a high 
density inner planet and a lower density outer planet, as 
is found in the Kepler 36 system. The outer planet could 
continue to accrete gas, increasing its mass and lowering its 
mean density after the swap took place. 

These integrations illustrate that the 'planet trap' set- 
ting, with numerous embryos embedded in the gas disk ex- 
terior to two planets at the edge of a clearing in the disk, 
is promising as it could account for both the proximity and 
high density contrast of the Kepler 36 planets. 

Transit timing observations can place constraints on or- 
bital inclinations (as discussed by Carter et al.|2012 1. When 
collisions between planets and embryos occur, planet incli- 
nations can be excited. For example, at the end of the inte- 
gration shown in Figure [3] the two planets differ by a degree 
in their orbital inclination. Future observations of the Kepler 
36 system may determine if there is a difference in inclina- 
tion between the two planets. The resonant chain shown in 
Figure [4] contains an embryo at a different inclination than 
the two planets. It is possible that such an object could exist 
in the Kepler 36 system but is not detected in transit. At the 
end of the integration in Figure [5] the two planets are at a 
different inclination than they were originally but we found 
that the two planets are are approximately coplanar (their 
mutual inclination is less than 0.2°), so not all integrations 
ended with mutually inclined planets. In this simulation no 
embryos were ejected so the variation in angular momen- 
tum is due to grazing encounters. Our integrations do not 
conserve angular momentum during collisions as we do not 
record or adjust planet rotation rates. 




50000 100000 150000 200000 250000 300000 

time (periods) 

Figure 3. Stochastic forcing from planetary embryos can knock 
two planets out of resonance and in some cases allows them to 
swap locations. The integration displayed is from the X-series 
with parameters shown in Table [2] The top panel shows semi- 
major axes, with error bars set according to the body's eccentric- 
ity so that they illustrate pericenters and epicenters of the orbit 
as a function of time. Here the red and green points are the two 
planets and the other color points correspond to planetary em- 
bryos. The middle panel shows the ratio of orbital periods for the 
two planets (outer divided by inner) and major mean motion res- 
onances arc shown with horizontal lines. The bottom panel shows 
inclinations of the bodies, with the same color points as the top 
panel. This integration illustrates that the positions of the two 
planets in some cases exchange or swap locations and leave the 
two planets in adjacent orbits. The planet originally closer to the 
star experienced no collisions with embryos. At the end of the 
simulation, the planet that experienced collisions is the innermost 
planet. If this planet was stripped during a collision, a dense core 
could have been left behind. This scenario could account for the 
large density contrast between the Kepler 36 planets. 



3.1 Properties of collisions 

For all collisions in both series of integrations we recorded 
the impact angles and relative velocities at impact of the two 
bodies involved in the collision. Histograms for the collision 
properties are shown in Figure|6]for planet/embryo impacts. 
The distributions are computed separately for each planet. 
The velocities are given in units of a particle in a circular 
orbit with the innermost planet's initial semi-major axis. 
The impact angle was computed with equation |23| We com- 



pute the critical grazing angle, (equation 24 1 using the radii 
from Table p] finding 9 g = 50° for impacts between a planet 
with the radius of Kep 36c and a planetary embryos with 
the same density but l/28th the mass, finding 9 g — 50°. 
Collisions with 9i m > 9 g are capable of stripping the man- 



tle of a planet ( Asphaug|2010 1 . The critical impact angle of 
9 g = 50° is shown as a wide, black, vertical line in Figure 
[6k. The impact angle distribution with both planets is wide, 



implying that both nornral and grazing impacts occur with 
both planets. 

As shown in FigurelSb, the impact velocity distribution 
is wider for the inner than outer planet. This is expected as 
embryos are likely to have higher eccentricities when they 
cross the orbit of the inner planet than the outer planet. 
We have shown the impact velocity distributions in units 
of the circular velocity. However the escape velocity from 
the planets is about 5 times larger than the circular velocity 
(see Table fT}. The mean of the impact velocity distributions 
are at a value of 0.2 times the circular velocity and so are 
approximately equal to the escape velocity of the planets. 
This value is in the regime (see Figure 8 by | Asphaug|2010 1 
where normal collisions are accretionary but grazing impacts 
are hit and run (little mass transfer or loss). The higher 
velocity impacts with the inner planet can cause disruption, 
and if they occur at grazing incidence then the envelope of 



© 0000 RAS, MNRAS 000, 000-000 




50000 100000 150000 
time (periods) 
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Figure 4. Similar to Figure [3] except illustrating a different 
simulation from the X-series. This integration illustrates that a 
resonant chain can be the end state. The middle panel shows 
that at the end of the simulation the two planets have the period 
ratio of the 7:5 resonance and the remaining embryo is in the 3:2 
resonance with the outer planet. 
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100000 150000 
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200000 



Figure 5. Similar to Figure [3] except showing an integration 
from the Y-series. This integration illustrates that the final result 
can be two planets in the 7:6 resonance. The innermost planet ex- 
perienced two collisions with embryos. If the planet was stripped 
during a collision, a dense core could have been left behind, ac- 
counting for the high density of Kepler 36b. 
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Figure 6. a) Distribution of collision impact angles measured 
from embryo/planet collisions in the X-series and Y-series inte- 
grations. Red and green histograms show collisions with the inner- 
most and outermost planets, respectively. Here an impact angle 
of 0° is a normal impact and an impact angle of 90° is grazing. 
Impact angles greater than 50° are sufficiently grazing that they 
can disrupt the envelope of a planet such as Kepler 36 c. b) The 
distribution of impact velocities in units of the initial innermost 
planet's circular velocity. Here a value of 0.2 is corresponds to 
the escape velocity of the Kepler 36 planets. The velocity distri- 
bution is wider for the inner planet. This is expected as embryos 
are likely to be at higher eccentricity when they cross the orbit 
of the inner planet. There is a wide distribution of impact angles 
and velocities implying that accretionary, disruptive, and enve- 
lope stripping collisions are possible. 



the planet can be disrupted leaving a dense core behind, as 
proposed for the origin of Mercury ( |Benz et al.|200f| . 

Figures [3] and [5] illustrate two possible scenarios for the 
formation of the Kepler 36 system. For the integration shown 
in Figure 13] the planet originally closer to the star experi- 
enced no collisions with embryos. All embryo/planet colli- 
sions occurred with the outer planet. However toward the 
middle of the integration, the two planets swapped location, 
and the planet that had experienced collisions with embryos 
became the innermost plant. If one or more of these col- 
lisions were sufficient high velocity and impact angle, the 
outer parts of this planet could have been stripped leaving 
behind the dense core that is currently Kepler 36b. Because 
the other planet was protected from collisions it would have 
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held onto its lighter elements and so would have remained 
at lower density. 

In the integration shown in Figure [5] that ended with 
two planets in the 7:6 resonance, the two planets did not 
exchange location. However, in this integration, the inner 
planet experienced two collisions with embryos. These col- 
lisions could have stripped the inner planet of low density 
material. Within the context of the planet trap, the outer 
planet could continue to accrete material. Its low density 
might in part be due to continued accretion. 

We did not allow mass to be stripped from planets dur- 
ing collisions in our integrations and did not take into ac- 
count orbital debris from collisions. The integrations that 
exchanged planet locations were from the X-series and so 
ended the integration with the more massive planet in an 
exterior orbit, opposite to the Kepler 36 system. Further- 
more in our integrations planets only gained mass from col- 
lisions with embryos. If an additional planetary embryo were 
to interact with the planets at the end of the integration 
shown in Figure [5] it is likely that the planets would not 
have remained in the 7:6 resonance, so the end states of 
the integrations could depend on embryo mass and number. 
Here we chose masses for the embryos large enough to strip 
a planetary envelope, and this was motivated by formation 
scenarios for Mercury ( Benz et al.||2007 |. However integra- 
tions with fewer more massive or more less massive embryos 
may also leave two planets in the 7:6 resonance. A larger 
number of more flexible integrations could be explored to 
determine which type of scenario would best account for the 
Kepler 36 system origin. 



4 DISCUSSION AND CONCLUSIONS 

In this study we have focused on origin scenarios for the 
Kepler 36 two planet system. This system contains two 
planets of very different density, in nearby orbits, that are 
just outside the 7:6 mean motion resonance. We first ex- 
plored a stochastic migration scenario for two planets of the 
same mass as the Kepler 36b, and c planets. Analytical es- 
timates show that precise adjustment of the migration rate, 
planet eccentricities or stochastic forcing parameter is re- 
quired so that two convergently migrating planets bypass 
lower j mean motion resonances, such as the 4:3 and 5:4, 
but capture into the 7:6 resonance. Both simple numerical 
integrations and a diffusive approximation for the stochastic 
forcing imply that once captured into the 7:6 resonance, the 
resonant system would not be long lived. Planets soon es- 
cape resonance and are likely to collide. Two planets could 
remain in the 7:6 resonance if the turbulent disk respon- 
sible for the stochastic forcing is depleted soon after reso- 
nance capture. Our findings conflict with the recent study 
by |Paardekooper et al.| ( |2013[ ) who have shown that long 
lived systems in the 7:6 resonance can be formed through 
hydrodynamic simulations that induce migration, eccentric- 
ity damping and turbulent stochastic forcing. Our analyt- 
ical estimate for the resonant lifetime and our numerical 
implementation of stochastic forcing are based on simplistic 
descriptions of stochastic diffusion so we may have under- 
estimated the resonant lifetime. Furthermore, the Kepler 36 
system is one of a few thousand planet candidate systems, 
and it is not necessary to avoid fine tuning in migration 



parameters or the disk depletion timescale as this type of 
system is not common. 

A mechanism involving stochastic forcing with disk tur- 
bulence alone would not account for the density difference 
between the Kepler 36 planets. We have explored a scenario 
motivated by the idea of a 'planet trap' (Masset et al. 2006 



Morbidclli et al. 20081. In this scenario two planets lie at 



the edge of a disk that contains a number of more rapidly 
migrating planetary embryos. We find that interactions be- 
tween planetary embryos and planets can nudge the two 
planets out of resonances such as the 4:3 and 5:4 resonances, 
leaving them in adjacent orbits. In our integrations with 
planetary embryos we have neglected stochastic forcing due 
to a turbulent disk, however interactions with a gas is still 
assumed to take place as we allow the embryos to migrate 
inwards. Integrations with a few approximately Mars mass 
embryos display a diversity of final states, including systems 
with two planets in the 7:6 resonance, resonant chains, colli- 
sions between both planets and planets that have exchanged 
location. 

We have recorded the properties of collisions that 
occur in the integrations. Approximately twice as many 
planet/embryo collisions occurred with the outer planet as 
with the inner planet and those with the inner planet had a 
wider distribution of impact velocities. For both planets the 
distribution of impact angles was wide, ranging from nor- 
mal to grazing impacts. The distribution of impact angles 
and velocities imply that the collisions can be accretionary, 
disruptive or strip the envelope of a planet. 

Two formation scenarios for the Kepler 36 planetary 
system are suggested by our integrations. An exterior planet 
that is stripped by a grazing collision with an embryo, could 
leave behind a dense core that then swaps locations with a 
protected inner planet, becoming a system with a dense in- 
ner planet in orbit interior to a low density outer planet, like 
the Kepler 36 system. Alternatively, an embryo could impact 
the inner planet and strip it in situ. The 'planet trap' dy- 
namical setting could account for both the proximity of the 
Kepler 36 planets as well as their large density difference. A 
different type of fine tuning is likely required as the num- 
ber and masses of interacting embryos are free parameters 
and the outcome is sensitive to the outcome of collisions. 



Motivated by a formation scenario for Mercury (Benz et 
ak]|2007 l, we chose masses for the embryos large enough to 
strip a planetary envelope. Simulations that include different 
numbers of embryos, a distribution of embryo masses, accre- 
tion and collisional disruption could be explored to better 
understand the origins of the Kepler 36 system. 

Our integrations show that when planets experience col- 
lisions with embryos, planet inclinations can be excited. As 
transit timing variations measure the masses of more planets 
it may become possible to search for correlations between 
planet densities and inclinations (and maybe even planet 
spin rates and obliquities) and so test the possibility that en- 
counters with embryos occurred in the late stages of planet 
formation. 
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